www.gusucode.com > 2D LINEAR TRUSS FINITE ELEMENT 工具箱matlab源码程序 > 2D LINEAR TRUSS FINITE ELEMENT/Linear_Truss_Finite_Element_2D.m

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 2D LINEAR TRUSS FINITE ELEMENT %%%
%%%     BY: CHRISTOPHER WONG.      %%%
%%%  UNIVERSITY OF NEW HAMPSHIRE   %%%
%%%      crswong888@gmail.com      %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% This program computes the displacements, support reactions, and stresses for a 2D truss finite element
%%% The formulation is for 2-noded elements with first order linear shape functions
%%% The SLE is solved for by eliminating i-th rows and i-the columns for which u(i)=0 (boundary conditions)
%%% A user is required to input truss information
%%% The algorithm should generally work for any truss structure


clear all  % clear workspace to avoid conflicts


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% PREPROCESSING (USER INPUT) %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%// Establish coordiantes and global dof's of nodes, node = [name, dof x, x, dof y, y]:
node = [ 1, 1,    0, 2,    0 ;
         2, 3, -500, 4,    0 ;
         3, 5,  400, 6, -300 ] ;  % x and y coordinates in mm
 
n_node = size(node,1) ;  % establish number of nodes (**NO INPUT**)
 
%// Discretization - assign elements to global nodes, elem = [name, node k, node l]:
elem = [ 1, node(2,:), node(1,:) ;
         2, node(1,:), node(3,:) ] ;
 
n_elem = size(elem,1) ;  % establish number of elements (**NO INPUT**)
 
%// Establish geometric and material properties:
E(1:6) = 70 ;  % Young's modulus (kN/mm^2)
A(1:6) = 200 ;  % Cross-sectional area (mm^2)
 
%// Establish applied loads:
P = sym('P',[2*n_node 1]) ;  % initiate load vector as symbolic array (**NO INPUT**)
P(1) = -20 ;  % load applied to dof 1 (kN)
P(2) = -19 ;  % load applied to dof 2 (kN)
 
%// Establish displacement boundary conditions:
u = sym('u',[2*n_node 1]) ;  % initiate displacement vector as symbolic array (**NO INPUT**)
u(3:6,1) = 0 ;  % dofs 3 through 6 pinned


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% GENERIC ALGORITHM FOR FEA OF 2D TRUSS %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%// Calculate length and the factors l and m for each element:
for i = 1:n_elem
    L(i) = sqrt( ( elem(i,9) - elem(i,4) )^2 + ( elem(i,11) - elem(i,6) )^2 ) ;  % length, mm
    l(i) = ( elem(i,9) - elem(i,4) ) / L(i) ;  % l factor (unitless)
    m(i) = ( elem(i,11) - elem(i,6) ) / L(i) ;  % m factor (unitless)
end
 
%// Construct global stiffness matrix from element local matrices:
K = zeros(2*n_node,2*n_node);  % initiate global stiffness matrix as 2D array of zero's
for i = 1:n_elem 
    k = E(i) * A(i) / L(i) * [     l(i)^2  l(i)*m(i)    -l(i)^2 -l(i)*m(i) ;
                                l(i)*m(i)     m(i)^2 -l(i)*m(i)    -m(i)^2 ;
                                  -l(i)^2 -l(i)*m(i)     l(i)^2  l(i)*m(i) ;
                               -l(i)*m(i)    -m(i)^2  l(i)*m(i)     m(i)^2 ] ;  % element local stiffness matrix (kN/mm)
                     
    index = [ elem(i,3) elem(i,5) elem(i,8) elem(i,10) ] ;  % global indices (dofs) for local stiffness matrix
    K(index,index) = K(index,index) + k ;  % add local stiffness matrices to respective global indices (kN/mm)
end
 
%// Establish which i-th rows and i-th columns will be kept for the subsequent calculation of 'u':
not_eliminated = false(2*n_node) ;  % initiate a false (0) logical array corresponding to the size of K
for i = 1:2*n_node
    if u(i) ~= 0  % condition for which a row or column is NOT eliminated, i.e., kept
        not_eliminated(i) = true ;  % set the i-th logical entry to true (1)
    end
end
 
%// Solve for unknown displacements by using the elimination approach with 'not_eliminated' rows and columns:
u(not_eliminated) = inv(K(not_eliminated,not_eliminated)) * P(not_eliminated) ;  % displacement solution (mm)
u = double(u) ;  % convert symbolic array to double precision
 
%// Solve for unknown reaction components of the vector 'P':
P = K * u ;  % includes applied loads and reactions (kN)
P = double(P) ;  % convert symbolic array to double precision
 
%// Determine elemental strains and stresses:
for i = 1:n_elem
    epsilon(i) = 1 / L(i) * [ -l(i) -m(i) l(i) m(i) ] * [  u(elem(i,3)) ;
                                                           u(elem(i,5)) ;
                                                           u(elem(i,8)) ;
                                                          u(elem(i,10)) ] ;  % element strain (unitless)
    sigma(i) = E(i) * epsilon(i) ;  % element stress (kN/mm^2)
end